Peripheral mixing of passive scalar at small Reynolds number 
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Mixing of a passive scalar in the peripheral region close to a wall is investigated by means of 
accurate direct numerical simulations of both a three-dimensional Couette channel flow at low 
Reynolds numbers and a two-dimensional synthetic flow. In both cases, the resulting phenomenology 
can be understood in terms of the theory recently developed by Lebedev and Turitsyn [Phys. Rev. 
E 69, 036301, 2004]. Our results prove the robustness of the identified mechanisms responsible for 
the persistency of scalar concentration close to the wall with important consequences in completely 
different fields ranging from microfluidic applications to environmental dispersion modeling. 



Mixing of passive tracer in a turbulent flow is a fun- 
damental problem of great technological interest which 
has undergone a significant theoretical progress in the 
last years IH, 0] ■ At low diffusivity and small scales the 
viscous-convective Batchelor regime of mixing arises for 
a large value of the Schmidt number Sc = v/k, where 
is the kinematic viscosity and k the molecular diffusivity 
p. In this regime the velocity field can be considered 
smooth because of the exponential decay of the veloc- 
ity spectrum. An analogous situation is realized if the 
tracer is advected by a time dependent, chaotic flow, i.e. 
a spatially smooth flow at all scales. Several theoretical 
predictions, including the exponential decay in time of 
passive scalar fluctuations 0, 0] , have been verified 
in experiments and numerical simulations Q . 

Recently, it has been shown that the presence of 
boundaries can alter significantly the predictions based 
on an unbounded domain. No slip boundary condition 
for the velocity field reduces the efficiency of mixing close 
to the boundary - the peripheral region - which becomes 
a source of passive scalar. This effect is of particular 
importance in the case of boundary dominated geome- 
tries, such as in microfluidics. Indeed, it is well known 
that the lack of efficient mixing is one of the problems 
in many microfluidic devices which operate at vanishing 
Reynolds number [l^]. Several statistical predictions for 
the peripheral mixing have been recently made [Til . Il^ 
and checked with laboratory experiments in a chaotic 
microchannel [l3l |. in polymer solutions [lA] and in kine- 



matic simulations 15 1 



In this letter, we will discuss the problem of tracer 
mixing in presence of boundaries by means of direct nu- 
merical simulation of a Couette flow at small Reynolds 
number and in kinematic simulations of a chaotic flow. 
Before discussing quantitative results, we can have a 
physical intuition from Figure [1] which shows a snapshot 
of tracer concentration advected by a two-dimensional 
chaotic flow. Because of the vanishing velocity, the tracer 
persists for very long time close to the boundary from 



FIG. 1: Snapshot of passive tracer concentration close to a 
wall placed on the top. White correspond to high concentra- 
tion, black to low concentration. The initial condition is a 
step function in the vertical direction. Tracers is advected ac- 
cording to IT]) with a two-dimensional synthetic velocity fleld. 



which it is only intermittently transported into the bulk 
in the form of elongated filaments. 

A passive scalar 9 advected by an incompressible ve- 
locity field u obeys the equation 



de 



u • V6' ^ kV' 



(1) 



where k is the molecular diffusivity, and appropriate ini- 
tial and boundary conditions are set for 9. In general u is 
subject to its own set of equations, the typical case being 
Navier-Stokes equations, together with boundary condi- 
tions and stirring mechanism which therefore determine 
the precise form of u close to the boundary. However 
some general predictions on the evolution of 9 can be 
made on the basis of the no-slip conditions and incom- 
pressibility only, assuming the flow to have a short cor- 
relation time with respect to the typical mixing time. To 
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simplify the notation, but without loss of generality, we 
assume the flow to be two-dimensional, with (x, y) the 
coordinates parallel and normal to the wall (y = corre- 
sponding to the wall) and (u, v) the associated velocity 
components. 

As a consequence of no slip (u(a;, y = 0) = 0) and 
incomprcssibility (V • u = 0) conditions, there is a 
region close to the wall, characterized by the scaling 
{v?'{y)) ^ and {v'^{y)) ~ y** (where by < • > we indi- 
cate the double average with respect to the coordinates 
parallel to the wall and velocity realizations). The ve- 
locities in the bulk of the container or duct are therefore 
much more intense with respect to the ones in the lay- 
ers close to the wall, and the passive scalar evolution 
becomes faster and faster as one moves away from the 
wall. Starting from these considerations, it is possible to 
describe the evolution of a passive tracer initially concen- 
trated in a layer of thickness 5 close to the wall in term 
of a turbulent diffusivity. Averaging ([T]), the equation for 
the y-evolution of the scalar profile is 



dt{e)^tidy[y%{e)] + ndl{e). 



(2) 



The first term in the rhs of ^ describes the role of 
chaotic advection in terms of an eddy diffusivity D = 
/o {"viv, 0)f (y, t))dt ~ /iy^. Comparison with the second 
term in the rhs of ([2]) suggests that the evolution of the 
profile is dominated by advection as long as (5 > ru = 
(k//^)^/'', the thickness of the diffusive boundary layer. 
Under this condition, the diffusive contribution can be 
neglected and ^ becomes 



dt{e) ^ tldyy^dyiO) 



(3) 



Taking as initial condition for the scalar a distribu- 
tion concentrated at the wall with 9{x, 0; 0) = 1 and 
limy^oo(^{x,y;t)= 0, the asymptotic solution of ([3]) for 
large times is [12 1 
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(4) 



i.e. the profile has a universal form, independently on 
the details of the initial distribution. The thickness 
S = {^t)^^/^ is the only characteristic scale and decreases 
in time, as the layer occupied by the scalar contracts in 
the evolution. It is important to remark that the specific 
form of the profile ^ gives a practically constant concen- 
tration for y < (5/4, making the boundary conditions for 
the scalar irrelevant in the advective stage. We remark 
that although ^ obviously conserves the average scalar 
(6), from Q one has that / {0{y, t))dy = 5(t)/^ ~ r^/^ 
is time dependent. The reason is that in deriving ^ the 
bulk is considered as an infinite reservoir for the scalar 
which has therefore zero average. 

Neglecting diffusion, the advection equation ([1]) holds 
for any local function of too and therefore its average 



is governed by a generalization of ([3]). In particular the 
moments of scalar concentration (0") are expected to fol- 
low the same profile (0]) independently of n. This is the 
mathematical consequence of the intermittent nature of 
scalar advection shown in Figure [1] in which (0") is dom- 
inated by the white regions in which 9 = 1. It is also 
possible to derive an expression for the probability den- 
sity function valid for any time and distance from the 
wall as 
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where g{x) is a Gaussian distribution with zero mean 
and variance 2/it and the dependence on 9 is defined im- 
plicitly through the monotonic, homogeneous along the 
wall, initial profile: 0o(yo) = & (note that large 9 corre- 
spond to small yo). At variance with which is valid 
only asymptotically in time, the prediction for the PDF 
is valid for any times as it depends on the details of the 
initial distribution. It is interesting to observe that ([5|) 
predicts that P{9 = 1, y, t) = for any y > 0: in the ab- 
sence of diffusion it takes an infinite time for the scalar 
to be transported away from the wall. 

We have tested the above theoretical predictions by 
means of numerical simulations of a three-dimensional 
Couette flow and a synthetic two-dimensional flow. For 
the first case, the Navier- Stokes equations have been nu- 
merically integrated in a 3D slab geometry of dimension 
Lx X 2Ly X Lz with no-slip boundary conditions on the 
two planes y = and y = 2Ly, and periodic boundary 
conditions on the streamwise and spanwise directions x 
and z. The flow is forced by the relative motion of the two 
opposed walls with opposite velocity ±/7o, from which a 
large scale Reynolds number is defined as Re = UoLy/v. 
Direct numerical simulations are performed by means 
of a standard pseudo-spectral Fourier-Chebyshev code 
Til . 17 1 at resolution 128 x 65 x 128 for a domain of size 



xL^ = 8x2x8. We use a moderate Reynolds 
number Re ~ 600, which is sufficiently large to sustain a 
turbulent-like motion for long time [l8'| but still small in 
order to have a well developed viscous layer where scal- 
ing imposed by boundaries is observed. The average wall 
shear rescaled with mean shear is so/{Uq/ Ly) ~ 3.3, the 
rescaled friction velocity u^/Ui^ = ^vsq/Uq ~ 0.074 and 
the friction Reynolds number Re* ~ 44. Scales and times 
are made dimensionless with the half-channel height Ly 
and large scale time T = Ly/Uo. 

In order to avoid the effects of diffusivity and increase 
the available observation range, the simulation of the 
tracer advection |T]) was carried out with a Lagrangian 
method: trajectories of = 10^ particles, representing 
the concentration of tracer, are integrated according to 
the equation x(t) = u(x, t). The initial condition for the 
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FIG. 2: Profile of tlie scalar density 9 at diflferent times in 
the Couette channel rescaled with respect to S{t) and com- 
pared with the theoretical prediction Q (continuous line). 
The inset shows the fitted values of 5 at different times and 
the prediction S ~ t~^^^. 
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FIG. 3: Average profile of the scalar density 9 at different 
times rescaled with respect to S and the bulk average b (see 
text) and compared with Q (continuous line). In the inset, 
the thickness of the profile is shown together with the predic- 
tion t~^^'^ dependence. 



particles is an uniform distribution in the x and z direc- 
tions in a layer close to the wall y < O.OALy. Continuous 
tracer distribution 0(x, t) is reconstructed at every time 
proportional to local tracer density on the grid. The ad- 
vantage of the present method is the possibility to per- 
form the simulation at Schmidt number virtually infinite 
(although a small numerical diffusion is always present), 
at the price of some noise in the reconstruction of small 
scales due to the discreteness of the tracer. 

In Fig. [2] we plot the mean profile {9{y,t)) at different 
times compared with the theoretical curve ([4|) . The value 
of the thickness 5 at different time is obtained from the 
fit of the profiles with ^ and its dependence on t is 
compatible with the prediction 6 = (/xi)^^/^, from which 
/i ~ 3.2 is estimated. 

Figure [2] shows that scalar profiles in the bulk deviate 
from the theoretical curve at large y. This is because, 
even if the value of Re is small, the scaling character- 
istic of the viscous layer, can be observed only up to a dis- 
tance y O.lLy from the wall. Therefore, in order to ex- 
tend the range of scaling, we performed an additional set 
of simulations based on a kinematic model for the velocity 
field. We define a two-dimensional velocity field in terms 
of a synthetic stream function \['(a;,y) = ^{x,y)G{y), 
where $ = sin(fca;a; -|- (px{t)) sm{kyy + (py{t)) represents 
a time-dependent cellular flow while G{y) is tailored to 
reproduce the correct scaling at the wall. If G ^ y^ for y 
close to the wall, the scaling of both components of the 
velocity is ensured. We chose G{y) ~ [l — cos{ksy))] close 
to each wall, with kg = 7r/(2Ls) and the position Ls of the 
inflection point defining the width of the scaling region 
for the velocity, y Lg, which is the region of interest, fn 
the bulk, a matching function connects the profiles. The 
phases of the cellular flow, Lpx and Lpy are given by a ran- 
dom process with a flnite correlation time. The velocity 



is similar to that used in 
approximatively to Ls = 



field generated by $ is placed on a grid of size = tt and 
2Ly = An at resolution 512 x 2048 where the evolution of 
^ is integrated by means of a pseudo-spectral code, with 
periodic boundary conditions. This numerical approach 

Ull^. Scaling regions extend 
4. The time unit is chosen so 
that the correlation time of the velocity field is T = I. 
In these units we have fi ~ 2.66 and k = 3.42 x I0~^ 
and therefore the width of the diffusive boundary layer 
is ru ~ O.G34Ly. As initial condition we choose a dis- 
tribution null in the bulk and concentrated at the walls 
in two smoothed-step functions of size Ly/A. The results 
are based on ensemble average over 100 realizations of 
the random noise driving the kinematic velocity field. 

In Fig. [3] the profile {6{y,t)) is compared, at different 
times, with the theoretical curve. In order to accurately 
resolve the region close to the wall, the extension of the 
domain in the bulk is not large enough for the approxima- 
tion of an infinite bulk to be valid. Therefore, because of 
the conservation of (9) , after a short transient a relevant 
amount of scalar accumulates in the central region of the 
domain, thus affecting the overall shape of the profile. In 
order to compare the numerical results with the theoret- 
ical prediction based on an infinite basin, we computed 
the profile of the auxiliary field 6 — {0 — b)/{l — b), where 
b{t) is the time-dependent value of 6 in the bulk (aver- 
aged over x). Figure [3] shows the remarkable agreement 
obtained between theory and numerics, indicating that 
the profile (jl]) can be easily extended to the general case 
of advection in a finite vessel. From the fitting proce- 
dure we get the values of the parameter 6, which is found 
to follow accurately the prediction S{t) = {fj,t)~^/'^ with 
/i ~ 2.13 (see inset of Fig. [3|). 

Figure [3] shows the profiles of different moments of 
scalar concentration (0"(j/, i)) computed at an interme- 



4 




FIG. 5: PDF's of the scalar density in the kinematic simula- 
tion at y — Ly/4:, corresponding to the thickness of the initial 

FIG. 4: Average profiles of moments of scalar density (61") distribution at times t = 0.02 (+) and t = 3.0 (x). Lines 

for n = 1 (+), n = 2 (x), n = 4 (*) and n = 6 (□) at time show theoretical predictions. 

t = 18T in the two-dimensional kinematic simulation. 



diatc time. All the moments collapse on the prediction 
(|4]), confirming the fact that in this stage diffusion is 
negligible and mixing of the scalar is dominated by eddy 
diffusivity according to ([3]). 

We have also computed the probability density func- 
tion, P{9,t,y), of the passive scalar. The results, at a 
distance y = Ly/A from the wall, corresponding to the 
thickness of the initial distribution, are shown in Fig. [5] 
for two different times. At very short time the PDF is 
peaked around the initial condition ??o/2. Observe that 
the distribution is much more narrow than that theo- 
retically predicted by (O, because the eddy diffusivity 
approximation ([2]) is not justified at short times. 

At longer time the PDF forms two pronounced peaks 
at the extreme values which are responsible for the sat- 
uration of the moments of 6, since the peak ui 6 — 1 
gives a dominant contribution to any (0"). Physically 
the presence of the two peaks is related to long tongues 
protruding from the wall region into the bulk. Advection 
stretches such structures, while preserving the value of 
0, the smoothing effect of molecular diffusivity becoming 
effective only on longer times. 

In conclusion, our results for both the three- 
dimensional channel fiow and the two-dimensional syn- 
thetic flow are accurately explained by the theoretical 
description of Ref. [l2|. Our results emphasize the im- 
portance of a correct description of the close-to-the-wall 
region where a scalar field tends to persist. The incor- 
rect reproduction of near-wall scaling behavior u{y) cx y 
and v{y) cx y^ necessarily destroys the above mechanism. 
This might be a serious problem in applications of en- 
vironmental dispersion modeling under strong stability 
conditions where the viscous layer can reach values of 
orders of meters thus affecting the realm of human activ- 
ities. A poor (in term of resolution) description of this 
region would result in a dangerous underestimation of 
the level of pollutants concentration close to the ground. 
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